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Abstract 

We consider active-sterile neutrino oscillations in the early universe in an inhomogeneous 
isocurvature background. We show that very small initial baryonic seed-inhomogeneities 
can trigger a growth of very large amplitude spatial fluctuations in lepton asymmetry. Do- 
mains of varying asymmetry are observed to persist for a long time despite dissipative effects. 
The space dependent asymmetry profiles give rise to MSW-resonances within the domain 
boundaries, enhancing dramatically the equilibration of the sterile neutrino species. Accord- 
ing to our one-dimensional toy-model, the effect is so strong that almost the entire parameter 
space where exponential growth of asymmetry occurs would be ruled out by nucleosynthesis. 



1 Introduction 



The present observational situation in neutrino physics appears to favour a conventional 
explanation of both solar neutrino problem and the atmospheric neutrino problem with large 
angle mixing among the ordinary, known neutrino species |l]. However, these scenarios leave 
unexplained the anomaly observed in Los Alamos neutrino experiment. To account for 
all known anomalies one necessarily must invoke new, sterile neutrino states u s , mixing with 
ordinary neutrinos. Such states also occur naturally in many particle physics theories beyond 
the minimal standard model. Active-sterile mixing would have several very interesting effects 
in astrophysical settings pfl and in particular for the evolution of the early universe 0- 



I7|. For example, the otherwise inert sterile states would interact with ordinary matter 



through mixing induced oscillation effects, and under the most natural assumptions could, 
for a wide range of oscillation parameters, be brought into thermal equilibrium prior to 
nucleosynthesis. The ensuing increase in the expansion rate of the universe could then bring 
the nucleosynthesis predictions for light element abundances in conflict with the observations 
in particular by causing an overabundance of helium-4. This would in particular 
be the case for the active-sterile neutrino oscillation solution for the atmospheric neutrino 
problem ||. 

Even outside the region of mixing parameters where equilibration is effective, one en- 
counters other interesting phenomena. For negative mass squared difference 5m 2 < 0, the 
lepton asymmetry L has been shown to evolve into an instability, which triggers an ex- 
ponential growth of L at the resonance temperature T res M. This idea has been used to 
circumvent the abovementioned exclusion of the — v s mixing solution for the atmospheric 
anomaly jlt|. The idea is that prior to temperature at which v s would be brought into 
equilibrium by z/ M — v s oscillations, a resonance in another mixing sector, say v T — v' s , creates 
a very large lepton asymmetry (without equilibrating u' s ), which then suppresses the am- 
plitude of oscillations in the — z/ s -sector enough to keep v s out of equilibrium until after 
nucleosynthesis. 

Later it was shown that the asymmetry growth after the resonance is chaotic in the sense 



that the final sign of the asymmetry cannot be simply deduced from initial conditions [11 
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It was recently verified |T^,|T^] that this effect is real, and that the sign is sensitive not only 
to the initial conditions, but for a wide region of oscillation parameters ( "chaotic"[] region) 
is also is unstable against small changes in parameters themselves^. Large asymmetries in 
electronic sector can significantly alter helium-4 abundance, which then cannot be reliably 
predicted should oscillation parameters ever be observed to lie in the chaotic region |T8|| . 

In this letter we further study the active-sterile neutrino oscillations and exponential 
growth of the asymmetry. Here we induce the additional, and completely natural, ingredient 
of spatially inhomogeneous initial conditions. Indeed, many baryogenesis mechanisms as 
well as possible first order QCD phase transition would predict isocurvature fluctuations in 
the baryon asymmetry in the early universe which persist all the way to the nucleosynthe- 
sis. What makes these fluctuations interesting for us is that they can provide seeds for a 
large inhomogeneous growth of lepton asymmetry, given the resonance conditions and the 
sensitivity of the direction of growth on initial conditions. We will show that this is indeed 
the case and that even very small baryonic seed-inhomogeneities may trigger a growth of 
very large amplitude fluctuations in lepton asymmetry. It should be noted however, that 
the occurrence of large amplitude fluctuations in L is not confined to the chaotic region of 
parameter space. We observe domain creation also in simulations involving "stable" mixing 
parameters, where the final sign is robust and defined by the sign of the effective asymmetry 
at T res . The domain structure in this case is triggered by the combination of inhomogeneous 
baryon asymmetry in the background and the diffusion effects, which together cause the 
effective total asymmetry seen by neutrinos to fluctuate in sign at the time of the resonance. 
This mechanism was first discussed in f22|| . 

The domains of varying asymmetry are seen to persist for a long time, the semistabil- 
ity being caused by a strong local MSW-effect which tends to restore large L against the 
smoothing effect of diffusion. Moreover, the MSW-resonances within the domain boundaries 
lead to a rapid equilibration of the sterile neutrinos. In our one-dimensional model this 



1 With the word chaotic we do not mean chaoticity in the mathematically strict sense, but that the 
solutions are highly sensitive to both initial and boundary conditions. 

2 There has been some debate as to whether the sign is chaotic or not. The analysis of ref. |2(J has recently 
been shown to be flawed |^,[^J] while the most recent numerical studies a momentum dependent quantum 
kinetic equations are not in conflict with the results of (l8| . 
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effect is so strong that the sterile neutrinos are brought close to equilibrium over almost the 
entire region of mixing parameters where the instability leading to the L-growth has been 



observed [12| 13 . Would this results hold also in the realistic 3-dimensional case, it would 



clearly invalidate the scenario of avoiding v s — equilibration, put forward in refs. [[10 
and discussed above. In short, the remedy would not work because in the first resonance 
involving z/ s /, along with building up the asymmetry needed to prevent ^-equilibration, the 
species v s < itself would be equilibrated. However, in exchange, one would obtain new, much 
stronger BBN-constraints on mixing parameters than one finds in the spatially homogeneous 
computation |J. 

We cannot however draw firm conclusions based on our 1-dimensional calculations. In 
fact, diffusion can be expected to be more efficient in higher dimensions, and also the realistic 
momentum dependent MSW-effect could be somewhat weaker in upholding the domains. 
However, it is clear that without a very careful account of the effects of inhomogeneities, the 
BBN constraints on mixing parameters may be underestimated, and moreover, the feasibility 



of the mechanism of |T(| remains in doubt. 

In section 2 we will set up our formalism by a variant of a moment expansion for the 
quantum kinetic equations valid for slowly varying semiclassical background. The outcome 
is a slight generalization of the single momentum approximation for QKE's with the in- 
clusion of diffusion corrections. In section 3 we detail our numerical approach, solve the 
evolution equations and present the results in our 1-dimensional model. Section 4 contains 
our conclusions and outlook. 

2 Diffusion equations 

Our starting point is the generalization of the usual Boltzmann equation for the particle 



distribution function to the case of a density matrix of a mixing 2-component system 



dtpij + ^{d P H, d*p}ij - ^{d x H, d p p}ij + i[H, p\ {j = C^p]. (1) 
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The Hamiltonian H of the system (we have dropped an irrelevant constant diagonal piece) 
is in the flavour basis given by 

H = —Ueo.Ul + V a a + , (2) 
2p 

where 5m 2 = m 2 U2 — , a + = (1 + cr 3 )/2 and # is the vacuum mixing angle. is the 
usual 2x2-vacuum mixing matrix with the convention that the mass eigenstate v\ becomes 
the active state in the limit 9 — > and V a is the matter induced effective potential for the 
active species. The curly brackets in ([I]) represent anticommutators. These terms are a 
generalization of the flow derivatives of a scalar Boltzmann equation to a case of a mixing 
multicomponent system, guided by the principle of hermiticity. The commutator term com- 
pletes the Liouville operator on the l.h.s. The collision term CV, [p] has the components (no 
sum over repeated indices is implied) 

cm = -d(i - <y^. - r el a±( PiJ - Pcq ) + cf\ P i (3) 

where we wrote the elastic collision term in the relaxation time approximation and separated 
from it the purely off-diagonal quantum damping matrix with D = \^ e i- The function p eq (p) 
is proportional to the equilibrium distribution for massless fermions, with the property that 
the integral over the elastic interaction term vanishes. Finally C inel [p] is the standard inelastic 
collision integral, and it is also proportional to a + . 

The set of quantum Boltzmann equations ([T]-|3]) provide an adequate description for a 
coupled mixing system of neutrinos in the presence of decohering scatterings and varying 
background, as long as the the length scale of the background variation is much larger 
than the compton wave length £c of particles. Here this condition is satisfied by an ample 
margin, since if, is some fraction of the Hubble radius whereas lc ~ 1/p ~ 1/T. 

We now make several simplifications. First, we observe that in the relativistic limit 
H ~ p, so that to a very good approximation Eq. (|]) reduces to 

d t pij + p • d x pij + i [H, p}ij = dj [p] , (4) 

where p = p/p. This equation is still a very complicated to solve, and we will proceed by 
reducing it to a truncated set of moment equations. To this end we make the following 
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decomposition 

pip) = PUM + 5p(p); I ^ pip) = 0, (5) 

where p = p{(p)) and the integral condition follows from a normalization assumption J p = 
pN eq . Integrating Eq. (|J) over momentum, weighted with zeroth and first power of velocity 
(here we have simply v = p = p/p), gives the following two equations 

dtPij + <9 X ■ Ay + i[H, p}ij = (dj) (6) 
<9t Ay- + d x pij + i[H, A]ij = (pdj), (7) 

where the thermal averages are defined as (• • •) = f ■ ■ • f eq (p)/N eq , and we used the shorthand 
notations for the averaged Hamiltonian H = (H) and the displacement vector A = (p5p). 
We also made the usual factorization approximation ([H,pSp]) — [H, A]. It is straightfor- 
ward to see that the moments of the collision integral appearing in (^|- 17|) are given by 

(Q 3 ) = -(l-5 tJ )Dp tJ +a±T^(nl q -n Ua n Pa ) (8) 

(pa,-) = -((7-J + i(i-^))r d Ay, (9) 

where n Va = p aa and n eq is the equilibrium number density of a species of massless fermions. 
Equations (|6|-||) form a closed set of 12 equations (24 including antiparticles), which are 
already much more feasible for numerical solution than the initial full QBE's. Solving them, 
one would obtain Pij{t), which represents the evolution of the total number of neutrinos, 
together with Ajj(i), which describes the total spatial flow of neutrinos. These equations 
are, however, still quite demanding to solve numerically. We shall therefore make a more 
constraining, diffusive assumption, according to which the Liouville operator acting on A in 
Eq. (^) can be neglected in comparison with the collision term. We then have 

■ A,,- ~ - J-ajfo, , (io) 

J- ij 

where we have defined r y - = (ay + |(1 — ^^Teq. Using ([K]) we can eliminate A from 
Eq. ([?D, yielding a single scalar (diffusion) equation: 

dtpij + i[H, p\ij + (1 - 5ij)Dpij = ^dlpij + afjC aa , (11) 

1 ij 



were we defined yet another shorthand notation C aa = T aa (nl q — n Ua rip a ). Apart from the 
diffusion term ~ d^pij, Eq. ( [TT]) corresponds to the usual set of "one- momentum" evolution 
equations for the reduced density matrix |J. The diffusion term describes the smearing of 
the local fluctuations, generated by the resonant amplification of the asymmetry, due to the 
random motion of neutrinos. 

It should be noted that the diffusion approximation, which is very well satisfied in normal 
systems, such as mixing fluids for example, actually breaks down for sterile neutrinos whose 
interaction rate Tf s is very small (if not zero as we have assumed so far). Nevertheless, 
isotropically generated disturbances average out even in a strictly collisionless case due to 
free streaming (a feature clearly missing in the case of mixing fluids). We will use a small 
but nonzero value for Tf s to regulate our expressions 3 (in practice we will take r„ to be 
about a tenth of r^ Q ). This approach mimics the effect of free streaming, except that it 
smooths out random noise that would be present in a realistic case. Moreover, an effective 
diffusion length is effected on the sterile species even when T ss = by the combined effect of 
scattering in the active sector and the mixing induced by the Liouville term in (^). This is 



presumably the reason why our results based on (|TT|) do not depend strongly on the actual 
(sufficiently small) value chosen for T ss . 

It is convenient to parameterize the reduced density matrices of the neutrino and anti- 
neutrino ensembles in terms of the Bloch-vector representation: 

p„(x) = ~ (P (x) + P(x) • a) , p p (x) = l - (p (x) + P(x) • a) , (12) 



where a are the Pauli matrices. The coupled equations with the diffusion terms from Eq. ( |lT|) 
in the case of v T — v s oscillations then read (other cases can easily be obtained by a simple 
redefinition 0): 

P = V + dlP + V_d 2 x P z + C TT (13) 
P = V + d 2 x P + V„d 2 x P z + C TT 
p = V x P - DP T 



3 It should be clear that this regulator leads to very little change in the sterile neutrino production rate; 
we of course do not introduce any inelastic interactions for the sterile states. 
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+ 



+ 



v-%p + v + dlP z +c TT ] z + v TS dl? T 

VxP- DP T 

v_d 2 P + v + dl? z +c TT ] z + v TS dlv T . 



where x = dx/dt, Pt = P x % + Pyf an d C TT = T TT (n 2 — n VT n DT ). The diffusion coefficients 



- — — 

2 I f^ 1 



(14) 



TT SS , 



pel _i_ "Tel 



should not be confused with damping coefficients D and D. The numerical values of the 
rates appearing above are |§ T TT = f rr ~ 0.32G 2 F T 5 and D ~ D = §r£. ~ 1AG%T 5 . 
The rotation vector V in the particle sector is 



V = V x x + (Vb + VL) z, 



(15) 



with the components 



5m 2 

W) 

5m 2 



sin 2^ 



cos 20- 17.8G F iV 7 -^r 

'2Ml 



2<P> 
y/2G F N^ L. 



(16) 

(17) 
(18) 



Here _/V 7 is the photon number density and the effective asymmetry L in the potential Vj, is 
given by 

(19) 



L = --L n + L Ue + L„„ + 2L, T (P) = L + 2L„ T (P) 

in the case of an electrically neutral plasma. Here L n is the neutron asymmetry and the 
term L Vt introduces the coupling between the particle and antiparticle sectors: 



l -(P Q + Pz )- l -(P G + P z 



(20) 



Finally, the rotation vector for anti-neutrinos is obtained by simply changing the sign of the 
asymmetry: V(L) = V(— L). 
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The neutrino and anti-neutrino ensembles are very strongly coupled in Eq. (|14]) through 
the effective potential term Vl(L). This makes their numerical solution very difficult. In 
particular, a technical problem is posed by the large cancellation in Eq. (^Of), where two 
0(l)-terms cancel up to 10 decimal places to produce a number of order O(L ) initially, 
with only slight improvement of matters before the very final stages of calculation when L 
becomes very large. This obviously gives rise to a very dangerous loss of accuracy in the 
numerical solution, which however can be easily avoided by defining new "large" and "small" 
variables: 

P± = P,±P,. (21) 
With these definitions the evolution equations become 

P + = V + cP x P++V_d?P+ + 2C TT (22) 

P = V + dlP +V_d 2 x P- 

P+ = V X P+ + V_d 2 x P + + V + d 2 x P+ + 2C TT 

P- = v x p; + v_d 2 x p +v + d x p- 



Pt = -V P±-V L P?-DP± + V TS d x P± 

P± = V Q P± + V L P^ - DP± -V x Pt + V TS d 2 x P±, 

Note in particular that the difference P _ is n °t affected by collisions. This occurs because in 
the momentum averaged approximation P is proportional to the total number of particles 
which can only be changed by annihilations, which affect particle and antiparticle systems 
in exactly the same way (T aa = r SQ ). 

It is now well known that in the homogeneous case the Eq. ( P^D can lead to a period 
of exponential growth of asymmetry []4, 12, 13| after the instability sets in at the resonant 
temperature || 

T res ~ 16.0 (|<W| cos2fl) 1/6 MeV. (23) 

This temperature sets the overall scale of importance for neutrino-oscillations in the early 
universe. The exponential growth of L is the key factor also for our present study, because 
it can lead to amplification of the spatial variations in leptonic asymmetry, seeded by the 
fluctuations in baryonic asymmetry contained in the term Lq in Eq. fll9|). These fluctuations 



are the most crucial feature of our model, and they naturally arise via several mechanism 
in yet earlier stages of the evolution of the universe, such as electroweak and QCD phase 



transitions 25 



To be concrete, we assume that the baryonic domain size is roughly one hundredth of 
the Hubble radius at the time of the QCD phase transition, ds — 0.01/if ()!;qcd)- Inserting 
the value of the Hubble expansion rate H(tQcr>) and accounting for the expansion of the 
universe, it then follows that at the initial temperature T < m M of our calculation the size 
of each distinct lepton domain is roughly given by 

d L ~ 0.0018—^- . (24) 

1 QCD J- 

Moreover, we have we have assumed that in each domain L Q (x) has a fixed value, which 
varies randomly in the range [0.5 x 10~ 10 , 10~ 10 ]. This is very conservative assumption, and 
fluctuations with orders of magnitude larger amplitudes can easily be created during the 



QCD transition ||25|| . The domains created in electroweak transition are typically smaller by 
a factor of thousand in size, and hence too small for us to resolve with our present numerical 
methods. It should be stressed that these parameters are just particular physically motivated 
choices; neither the actual domain size or the range of amplitude variation have particular 
importance for our final results. All that matters is that some fluctuations exist to launch 
the instability. 

Solving equations (E3T) is numerically relatively easy in comparison with the full mo- 
mentum dependent equations, resulting in huge drop in the computer power required. We 
believe that they nevertheless provide a qualitatively correct description of the essential fea- 
tures of asymmetry oscillations in the spatially varying background. In fact, the numerical 
difficulties inherent in the full momentum dependent equations are so severe that at present 
the asymmetry oscillations can not be distinguished with certainty from computer-generated 
numerical errors even in the spatially homogeneous case | |26|| , and hence using them in the 
spatially varying case would appear an insurmountable task. 

Let us finally comment on feasibility for the present problem of the so-called static ap- 
proximation [|1^,|13|, which is often used to simplify the full momentum dependent quantum 
kinetic equations. The static approximation does include momentum dependence, but has 
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the great disadvantage that it averages out the asymmetry oscillations, whereby the effect of 
creation of leptonic domains in the chaotic region is poorly described in that approximation. 
More importantly, the static approximation excludes the MSW-resonance occurring in the 
boundaries of the domains. Since these resonances are responsible for the abundant pro- 
duction of sterile neutrinos we are finding in our computations, we conclude that the static 
approximation fails to describe the most striking feature of the asymmetry oscillations in 
the inhomogeneous backgrounds. Hence the only reasonable improvement over the present 
approximation is to go directly to full QKE's, which, as stated above, appears an overly 
difficult task to be implemented numerically. 

The main advantage of evolution equations (|23|), over the previous treatments [p2| , |2"T|j is 



the inclusion of the diffusion terms while keeping a full account of the asymmetry oscillations 
and the MSW-effect. This, we shall see, leads to new and interesting features in the solutions. 

3 Numerical results 

It is useful first to consider the information provided by earlier studies without spatial fluc- 
tuations as a function of the mixing parameters. Indeed, the (5m 2 , sin 2 2#o)-plane can be 
divided to regions with phenomenologically distinct characteristics as is roughly outlined in 
Fig. [I]. In the region left from the thick solid line there is no growth of asymmetry |T2"|,|T3 



This occurs because for very small mixing angles the resonance becomes so nonadiabatic that 
the solution has no time to pull away from the weak local fixed point at L — 0, which still 

one can 



exists below the resonance temperature for very small mixing. Using results of fT2 
estimate that this occurs for sin 2 2#o 15m 2 ! 1 / 6 < 10~ 10 . The asymmetry growth is thwarted 
also in the region to the right from the thick dashed line, because there the sterile states are 
brought to a full thermal equilibrium before the resonance M . A large asymmetry growth is 
observed in the rest of the parameter space. This region is further subdivided to three areas: 
the area above the thin solid line is called "stable region", since there the final sign of the 
neutrino asymmetry is determined by the initial baryon asymmetry. Below the thin solid line 
lies "chaotic region" , where the final sign of the neutrino asymmetry is very sensitive to tiny 
changes in either initial conditions, or oscillation parameters, due to a period of rapid sign 
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changing oscillations of the neutrino asymmetry right after the resonance. Slightly different 
results to the extent of the chaotic region exist in the literature: the region bounded by the 
thin solid line corresponds to the results of , whereas [ 26j reports a possibility of chaotic 
behaviour in a somewhat smaller area, extending to the right from the thin dashed line. This 
difference is likely attributable to the fact that [^] used the momentum dependent QKE's, 



whereas a momentum averaged equations were employed in fL8| . Finally, in the lower right 
corner below the dash-dotted line the asymmetry oscillations continue until the active neu- 
trino decoupling. Since our diffusive approximation is not valid anymore on this area, we 
shall not consider these parameters further in this paper. 

i 1 i 1 ^ 1 i 1 i 1 
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Figure 1: The mixing parameter space divided into regions with qualitatively different phe- 
nomenological aspects. 



Numerical setup. We have solved the set of partial differential equations Eq. (|23|) by the 
method of lines using central finite differences. Obviously the spatial discretization must be 
such that our grid encompasses all the relevant physical scales. On the small distance side 
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this implies that the grid must be dense enough to accurately model the initial domains of 
baryon asymmetry. Our numerical checks indicate that each domain should be divided to at 
least twenty sub-zones, to get sufficient accuracy in evaluation of the numerical derivatives 
in diffusion terms. On the other hand, the grid must be large enough to accommodate 
all diffusion scales (active and sterile) during the entire computation. This is particularly 
complicated, because the diffusion length increases as Id ~ T~ 5 as the temperature decreases. 
These constraints on the discretization lead to the need for very large lattices, which require 
a lot of computer power. Thus, for a first study, we have only considered one dimensional 
space. This is of course unphysical and to a degree compromises our ability to draw definitive 
quantitative conclusions, but we believe that our results provide important qualitative insight 
to the physical phenomenon under study. We will return to this issue in our conclusions. 

The complete solution of the problem requires setting the boundary conditions at the 
ends of the lattice. We have used periodic boundary conditions, which are very natural, 
since they automatically impose the desired conservation of particle number. One might 
think however, that periodicity is too restrictive when assumed also for the off-diagonals. 
We have hence investigated the effect of using various other possible boundary conditions 
on the system (cf. Fig. 3), but found no difference in results, which could not be attributed 
to finite size effects, as long as the particle number was conserved. 

To cope with the requirements on the grid size, introduced by the rapid growth of ip, we 
invented a technical trick of doubling the grid size periodically. The doubling step consists 
of simply discarding every second point in the original grid and replacing it with a double 
image of the resulting more sparse grid. This procedure can be repeated as many times as 
necessary, provided that no relevant physics is lost when coarse graining steps are taken. 
Certainly the initial baryon number fluctuations, whose scale remains frozen, are eventually 
evened out by the coarse graining. However, as will be explained later, these features are not 
relevant for the dynamics of the system immediately after the resonance. Indeed, we have 
checked the validity of the doubling method by comparing results of the codes running with 
and without a doubling step over an interval of time where both are expected to be valid, 
finding an excellent agreement. 
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We now describe the characteristic features of our numerical solutions. As in the spatially 
homogeneous case, the evolution of the asymmetry can be divided to three different phases, 
albeit the physics at each stage is markedly different. 



oscillations. This behaviour is very similar to the homogeneous case, except that the frozen- 
in baryonic fluctuations in the background asymmetry Lq now induce compensating spatial 
variations in the dynamical variable L Vt (P(x)). On the other hand, diffusion tends to smooth 
out these variations, thereby opposing the tendency of relaxation to £(x) = 0. The outcome 
is that just before the resonance the effective asymmetry L(x) will have spatial fluctuations, 
whose amplitude and scale are set by the initial amplitude and the size of the domains and the 
strength of the diffusion at the resonance. These fluctuations also typically oscillate in sign 
even though the initial fluctuations were all positive, as was recently noted by DiBari ||22|| . 
The sign fluctuations are requisite for the creation of the large amplitude fluctuations in L 
in the stable region of parameters. In the chaotic region the asymmetry generation is not 
dependent on this effect however, and there a spatially varying L Ut (P(x)) at the onset of the 
resonance is sufficient. In Fig. 2 we have plotted the variable L Ut (P(x.)) as a function of x 
for a small part of our grid in various temperatures for sin 2 29q = 10~ 5 and 5m 2 = — 1 eV 2 . 
The uppermost panel corresponds to a situation a little above the resonance, which for these 
parameters occurs at T res — 16 MeV. The distribution in this panel is effectively the negative 
of our initial baryon asymmetry fluctuation spectrum, smoothed out by diffusion. 

2. Resonance region. At the resonance the lepton asymmetry enters the region of in- 
stability. In other words, the previously stable fixed point L = becomes unstable, and 
asymmetry undergoes a period of exponential growth and sometimes violent oscillations 
(chaotic region). In the stable region, the domains of opposite effective asymmetry L(x.) 
start to grow rapidly in opposite directions. Similarly in the chaotic region, the spatially 
varying L Ut (P(x)) triggers a period of growth and incoherent oscillation of asymmetry in 
neighbouring regions. In both cases this behaviour is regulated by the diffusion, which tries 
to smooth out the fluctuations just generated. The oscillation induced asymmetry growth is 
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much stronger however, and asymmetry growth cannot be stopped. Second panel in Fig. 2 
shows how the spatial variations in L looks during the rapid oscillation period following the 
resonance this point. It should be noted that soon after the initial exponential growth period 
the dynamics of the system is solely determined by balancing diffusion against the MSW- 
effect, which is driven by the instantaneous asymmetry configuration; the initial asymmetry 
Lq plays no role at all, apart from triggering the initial growth period. 

-10 

1x10 

-6 

6x10 

-6 

6x10 

2x1 0" 5 
4x1 0" 5 

0.01 0.02 0.03 0.04 0.05 0.06 

Figure 2: Shown is the spatial variation of the neutrino asymmetry for five representative 
temperatures, indicated by the numbers to the left from the panels, for (sin 2 29 , 5m 2 ) = 
(1(T 5 , -1 eV 2 ). The units are £i 8 (18MeV/T) where £ 18 is the Hubble radius at T = 18 MeV, 
and the numbers to the right give the extent of the |/-axis in each subpanel, measured from 
bottom to the dotted line indicating L Vr = 0. Only a small part of the actual grid is shown. 




3. Annihilation region. After the resonance, once the asymmetry oscillations end and 
the asymmetry growth gets less violent in the chaotic region, diffusion effects take partly 
over. At this stage a quasistable distribution of leptonic domains with varying sign of L is 
created. The smallest domain size at this point is set by the diffusion length £jj of the active 
species, while a typical domain size is a few times diffusion length and some of the largest 
may be even ten times the diffusion length in size. The absolute value of the asymmetry 
in the domain centers approaches the new power law fixed point value of the homogeneous 
case : 

L fp oc ±T~\ T < T res , (25) 
14 



with |Lf p | ^> |Lo|. After this point the number of leptonic domains remains unchanged for 
some time, while the size of the domains adjusts slowly. Here one observes the peculiar 
phenomenon that the smallest domains do not vanish, but instead grow with the diffusion 
length at the expense of larger ones. This is caused by the extremely strong tendency of 
a local asymmetry to roll towards the fixed point value (|25|) . This "pull" gets stronger 
the further away from the fixed point L is driven, and thus, it is the relative size of the 
large domains that gives away first to diffusion rather than annihilating away the smaller 
ones. This stage of evolution is depicted in the third panel in Fig. ^. In the stable case 
large leptonic asymmetries are first created in boundaries domains with different baryonic 
asymmetry. The number and size of created domains is then initially sensitive to the size 
of the baryonic domains. However, when the diffusion effects kick in, the later evolution of 
the system is very similar to the one observed in the chaotic region; only the mechanism of 
domain generation is different. 

At some point, the smaller domains have eaten the larger domains resulting in a distri- 
bution of equal sized domains (fourth panel in Fig. 2), after which domains eventually begin 
to annihilate. However, as was explained above, there is a resisting tendency against domain 
annihilation, so that (in a finite grid) the system "supercools" against domain annihilation. 
As a result, when an annihilation eventually takes place, several domains are often destroyed 
at once (fifth panel in Fig. 2). In a finite grid such a multiple annihilation event shows up as 
a slight discontinuity in the sterile neutrino production rate. After an annihilation event the 
equalizing effect takes over and smooths the distribution again, until an even distribution of 
yet larger domains is again established, and the cycle takes over. It should be obvious that 
eventually the diffusion effect must win, completely smoothing out the asymmetry. However, 
our approximations break down well before this would happen and hence we cannot see this 
limit in our numerical examples. 

Physical consequences. The most important physical implication of the existence of the 
quasistable domains of varying asymmetry is that it gives rise to copious production of sterile 
neutrinos, as a result of an MSW-effect on the domain boundaries, induced by the spatially 
varying effective potential Vl oc L. Indeed, the large mixing at the resonance, combined 
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with the decohering interactions brings sterile neutrinos into equilibrium very efficiently. 
This mechanism was first discussed heuristically by Shi and Fuller [^7j. Here the mechanism 
is automatically built into our evolution equations. In Fig. [3] we show the ratio of the 
(sum of neutrinos and antineutrinos) active and the sterile number densities as a function of 
temperature for two representative choices of parameters. The upper figure corresponds to 
(sin 2 29q, 5m 2 ) = (10~ 8 , — 10 2 eV 2 ), taken from the stable region of parameters (the point A in 
Fig. 1 ), and the lower figure corresponds to (sin 2 2# , 5m 2 ) = (10 -5 , —1 eV 2 ) from within the 
chaotic region (point B). The rapid initial growth observed in each of the curves corresponds 
to the pre-inst ability near resonance when L pa throughout the space. Immediately after 
the resonance L grows rapidly and the MSW-effect becomes confined to within narrow bands 
inside the domain boundaries. Moreover, it occurs only for particles (L > 0) or antiparticles 
(L < 0) at a time. As a result the equilibration proceeds slower after resonance. 

The dashed and dash-dotted lines in the upper figure correspond to using different bound- 
ary conditions. For dashed line we imposed periodicity for the diagonals of p, but set the 
off-diagonals to zero at the boundaries. For dash-dotted line we imposed periodicity for all 
first derivatives of p, but did not restrict p itself. The difference in results is essentially a 
measure of the finite size effects, and become nonneglible only when the structure of equally 
sized domains is settling in. Slight "dents" seen in the upper figure, caused by the multiple 
domain annihilations, are another reflection of finite size of the grid. We have checked that 
these features become less prominent when the lattice size is increased, without changing 
our conclusions. We have also checked that the production of sterile neutrinos is robust 
to variations in the cut-off value of sterile neutrino interaction rate. In fact the effect of 
decreasing sterile neutrino interaction rate is a small increase in the rate of sterile neutrino 
production, which helps to justify our diffusive approximation. 

In the one dimensional case studied here the equilibration effect is very strong: while 
we have not tried to completely map the parameter space, it appears clear that for most 
of the parameter space where exponential asymmetry growth takes place, the sterile neu- 
trinos are excited with large enough energy density to be in conflict with nucleosynthesis 



constraints |[28|| . We have shown this explicitly in Fig. 3 for parameters corresponding to 
points A and B in Fig. 1. However, one does expect that for very small |5m 2 | no equili- 
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Figure 3: Ratio of sterile and active neutrino number densities as a function of temperature. 
Upper (lower) figure corresponds to parameters indicated by points A (B) in Fig. 1. The 
dashed and dash-dotted lines correspond to set A with different boundary conditions (see 
text). 

bration should take place, because then T res becomes very low, and consequently the overall 
scattering rate T ~ T 5 becomes very weak. We have checked that for sets C\ and C2 we still 
find relatively strong equilibration, whereas for sets D\ and D2, corresponding to T res ~ 5 
MeV, practically no //^-excitation takes place. So, below a line situated somewhere between 
these sets exists a region where a large L can be created without bringing v s into equilibrium. 



However, this region is already disfavoured by other reasons ||14|| . We have also checked for 
several parameter sets that for very small mixing angles (region left from thick solid line in 
Fig. 1) no sterile neutrino equilibration takes place. 

Another interesting result following from our studies is that, despite the large fluctuations 
in leptonic asymmetry, the average total asymmetry is found to be zero in our calculations, 
up to finite size effects. This is of course not surprising, but rather something one should 
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naturally expect even without a detailed computation; an average of a large number of sub- 
horizon scale fluctuations with random signs and symmetric amplitudes should be close to 
zero. Notably, such a configuration has much less effect for nucleosynthesis than would have 
a large scale homogeneous lepton asymmetry . 

Should our results hold also in a realistic 3-dimensional case, BBN would exclude the 



active-sterile neutrino mixing essentially in the entire interesting region of parameters [14 



However, one should expect that both the inclusion of the realistic momentum distributions, 
and the higher dimensional spatial structure could decrease the strength of our effect. Intro- 
ducing the momentum distribution would spread out the MSW-resonance condition on the 
wall boundaries, perhaps making it less efficient in opposing the diffusion effects. This might 
lead to somewhat shallower boundaries and perhaps easier domain annihilation, to which 
we already alluded above. In a higher dimensional lattice the diffusion would be much more 
efficient simply because of the topology; now the domains could spread out in more ways 
than one, and the tendency to create a quasistable domain configuration would be weaker. 
Unfortunately such simulations turn out to be too demanding numerically at the moment, 
and we cannot provide any quantitative estimate for the size of these effects. We therefore 
will only conclude that we have found, and verified in a toy model, an interesting new feature 
of asymmetry oscillations in the early universe, which has the potential to significantly affect 
(strengthen) the bounds on neutrino mixing parameters arising from BBN-considerations. 



4 Conclusions 

We have shown that large amplitude fluctuations in the lepton asymmetry can naturally be 
created in the early Universe. The main ingredients needed are the active-sterile neutrino 
mixing, which for a large range of mixing parameters encounters an instability leading to a 
period of exponential growth of asymmetry at T res ~ 0(10) MeV, and small baryonic seed 
asymmetries, which are naturally created for example during the electroweak and QCD phase 
transitions. The effect is not dependent on the particular values of amplitudes and distance 
scales of the seed asymmetries, which are only needed to trigger the initial inhomogeneous 
growth, because right after the initial growth period, the dynamics of the system is solely 



18 



governed by the MSW-effect and diffusion. 

The phenomenology of the system is completely different from what one finds in calcu- 
lations assuming a homogeneous background. First, the asymmetry does not reach a large 
homogeneous value during the BBN. Instead, a quasistable structure of sub-horizon domains 
with varying sign of asymmetry is realized, while the average asymmetry over the horizon 
scale remains close to zero. Such a configuration has much less impact on the light element 



abundances |1| , or cosmic microwave background |29| than a large homogeneous asymmetry 
would have had. Moreover, it is important to note that the creation of sub-horizon domains 
is a robust feature. In order to realize a scenario where only superhorizon scale perturbations 



are amplified [p^j , one must make the unnatural ad hoc assumption that the early universe 
is almost exactly homogeneous at all sub-horizon scales. If this were not the case, the tiny 
sub-horizon scale fluctuations would take over, giving rise to the structures discussed in this 
paper and thereby completely drowning any signature at superhorizon scales. 

An MSW-resonance would occur within the domain boundaries, where Vi oc L(x), in 
a manner very much similar to the inside of the Sun. As a result the mixing angle would 
be enhanced within a confined region inside each domain, which, together with the rapid 
decohering scatterings, leads to an efficient equilibration of the sterile species. We found 
this effect to be strong enough to lead into conflict with BBN-constraints over most of the 
region where the rapid asymmetry growth has been detected in homogeneous calculations. 
However, we stress that our computations are based on a one-dimensional toy model, and 
furthermore employ a momentum averaged approximation for the full QKE's. Relaxing 
either of these approximations might be expected to weaken the equilibration effect, although 
it is impossible to quantify how much. We conclude that the presence of baryonic seed 
inhomogeneities may have a major impact on sterile-active neutrino mixing in the early 
universe, and in particular on the BBN-bounds on mixing parameters. 
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